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Abstract. We perform a comparative analysis of different computational approaches 
employed to explore the electronic structure of ultralong-range Rydberg molecules. 
Employing the Fermi pseudopotential approach, where the interaction is approximated 
by an s-wave bare delta function potential, one encounters a non-convergent behavior 
in basis set diagonalization. Nevertheless, the energy shifts within the first order 
perturbation theory coincide with those obtained by an alternative approach relying 
on Green’s function calculation with the quantum defect theory. A pseudopotential 
that yields exactly the results obtained with the quantum defect theory, i.e. beyond 
first order perturbation theory, is the regularized delta function potential. The origin 
of the discrepancies between the different approaches is analytically motivated. 
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1. Introduction 

Diatomic ultralong-range Rydberg molecules consisting of a Rydberg atom whose 
electron, upon frequent scattering off a ground state atom, binds the atom, localized at 
large distances (~ 10 3 Bohr radii), to the ion core of the Rydberg atom were predicted 
theoretically in Ref. [T|. A subset of these Rydberg molecules with s-wave dominated 
electronic orbitals were later observed [2]. This initial experimental work was followed 
by a flurry of new observations of Rydberg molecules with a variety of electronic and 
rovibrational structures IMo). For electronic s-wave scattering of the Rydberg electron 
from the ground state atom, two types of ultralong-range molecular states can be 
distinguished: the non-polar, low angular momentum quantum defect states and the 
polar, high angular momentum “trilobite” states. Theoretically, some aspects of these 
high angular momentum states have been explored, notably the semiclassical nature of 
the rich oscillatory structure in the adiabatic potential energy surfaces m, their large 
permanent electric dipole moments Da. as well as the precise control over their electronic 
properties and molecular orientation by static electric and magnetic fields fl3]IT4j. 

Further emphasis has been laid on polyatomic systems consisting of a Rydberg atom 
bound to two or more ground state atoms dans] or to a diatomic polar perturber mm- 
Theoretical approaches to describe ultralong-range Rydberg molecules can be divided 
into two categories: methods using the Fermi pseudopotential mm and methods 
using the quantum defect approach [20]- for a recent overview, we refer the reader 
to ED- Although the two approaches are expected to be equivalent, they indeed differ 
in numerical implementation (T9] and are even used alternatively [3]. 

In this work we therefore aim at a concise analysis of the origin of the discrepancies 
obtained from these two approaches. We thereby point out the interconnection and 
limitations in both approaches which should be taken into consideration in future studies 
of ultralong-range Rydberg molecules. Specifically, we show that a pseudopotential 
modeled by a bare delta function potential leads to non-converging molecular potential 
energy curves in basis set diagonalization of the electronic Hamiltonian. Modeling the 
pseudopotential instead by a regularized delta function reproduces potential energy 
curves which agree exactly with those obtained by the quantum defect theory approach. 
Throughout our analysis we not only give limits of the bare delta function potential but 
also its validity as an approximation in first order perturbation theory. This ultimately 
links the two approaches. 

2. Ultralong-range Rydberg molecules 

We consider a Rydberg atom whose ionic core is located at the origin and a neutral 
ground state atom located at the position R within the Rydberg electron orbit. In the 
Born-Oppenheimer (BO) approximation, the electronic Hamiltonian for this Rydberg 
molecule reads 

H = H° + V(R) 


( 1 ) 
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where H° is the Hamiltonian describing the Rydberg electron in its ionic core potential 
and V (R) is the interaction between the Rydberg electron and the neutral ground state 
atom which we call from now on the perturber. If V is short-ranged, the corresponding 
electronic wave function T(r) outside the range of V can be determined from the 
scattering phase shifts induced by V. Although a more general treatment is possible, 
here we focus exclusively on pure s-wave electron-perturber scattering. 

In the region where the electron is near the perturber, but outside of the range of 
V, the electron wave function behaves as 


T(r) oc 


sin (p — 5(k)) 
kp 


( 2 ) 


where p — |r — R| is the relative electron-perturber coordinate and S(k) is the s-wave 
phase shift depending on the wave vector k given semiclassically as 


*+£)• < 3 > 
The s-wave scattering length of the potential V is obtained from a[k\ = — tan (S(k))/k, 
which in the low energy limit behaves as 

a[k\ = a[0] + —a[0 \k + 0(k 2 ) (4) 

o 

where a[0] is the zero-energy scattering length and a[0] is the zero-frequency (static) 
polarizability of the ground state atom (polarized in the charge-atom interaction) |2ll22] . 
The momentum k can be considered to be small (low-energy limit) because R in (J3]) 
is typically large and close to the classical turning points. The idea behind the Fermi 
pseudopotential approach is to replace V by a zero-range pseudopotential that possesses 
the same s-wave scattering length as the original potential, while the quantum defect 
approach replaces the potential V by the Dirichlet boundary condition (J2]). 


k = k(R, E) = J 2 



3. Bare delta function potential 

The pseudopotential usually employed to describe the s-wave binding of ultralong-range 
Rydberg molecules [1] is the bare delta function potential 

V (r, R) = 2Tra[k\5(r — R), (5) 

where the scattering length a[k] depends on the energy E and the perturber position 
R via (J3J) and (j4j). In the case that one is interested only in energies E close to some 
reference energy E® one may approximate 

k « k(R, £°) . (6) 

We now consider the energy shift of the electronic Rydberg orbit by the perturbing 
potential V where | is an eigenstate of the bare Rydberg Hamiltonian with eigen 
energy E°: 


iP\<p i ) = EP i \< Pi ) . 


( 7 ) 
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For a compact notation we choose the basis states | ipf) such that they diagonalize V 
inside energetically degenerate submanifolds, i.e. (ipf V \(pj) oc S i3 if E® = E { -. Using 
(l6j) . we deduce from (J5J) the hrst order perturbation theory correction to Ef as 

£f>(R) = 2™ [*(«,£“)] MR)| 2 , (8) 

which is the BO potential, valid for low angular momentum quantum defect, as well as 
for high angular momentum “trilobite” states. 

As an example we consider the Rb(35s ) + Rb(5s ) state studied in |2|,[3i|. The 
correspondent eigenstates of H° are the quantum defect states | nlm) that possess the 
energies E® lrn = — l/2(n —A;) 2 , where A i denotes the /-dependent quantum defect. We 
use A 0 = 3.13, Ai = 2.65, A 2 = 1.35 and neglect the quantum defects of higher, i.e. 
A;>2 = 0, angular momentum states. Here n, l and m are the usual hydrogenic quantum 
numbers. The properties of the perturber Rb(5s) enter in (|4j) via the dominating triplet 
s-wave scattering parameters 

a ra>[0] = —16.1 a 0 and «_Rb[0] = 319.2 Oq , (9) 

where a 0 is the Bohr radius. The e”-Rb(5s) triplet scattering length has been calculated 
in [23j and the polarizability is obtained from m 



Figure 1 . Comparison of the Born-Oppenheimer potential curves E as a function of 
the nuclear separation R for a Rb(35s) + Rb(5s) molecule with the parameters (O for 
different methods. The dissociation threshold is set to zero. E d) and E^f are the first 
order approximations based on (|8|). where uses the effective zero-energy scattering 
length a e ff[0] = — 18.5ao. The curve E leg (R) gives the energy for the regularized delta 
function potential based on (|20l) while the three curves Ej are determined via numerical 
diagonalization of CD including basis sets 32 — j < n* <32 + j. 
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The BO potential energy curve E d) is shown (red dashed line) in Fig. [TJ At 
R « 1900ao, it possesses a minimum that lies 26 MHz below the dissociation limit 
(here set to zero) and supports several bound states. The non-differentiability around 
R rs 2050ao occurs due to the non-analytic behavior of k and E at the classical turning 
points. Additionally, the dashed orange curve Fig. [[] shows the lower-lying potential 
curve based on the effective scattering length a e ff[0] = —18.5ao which was found 
out to agree better with the experimentally observed spectra [2]. As we will outline in 
Sec. [2 one reason for the discrepancy between a e ff[0] and a[0], which has not been fully 
understood so far, is the effect of couplings between non-degenerate \pf) and \tpf) basis 
states that, at the level of first order perturbation theory, are neglected in (J8]) but lead 
to deeper binding energies. This coupling becomes especially important in cases where 
the quantum defect of the state in question is close to an integer value. 

To take the couplings into account, which means to analyze the system beyond the 
first order perturbation theory, we diagonalize H numerically in an extended basis set. 
Technically this is done by fixing k as in (jBJ) to some energy E { - and diagonalizing El in a 
basis set of states range-bound to a certain number of quantum defect eigenstates | nlm) 
whose energies lie close to E®. In the case of Rb(35s) + Rb(5s), we choose E® = E® 5 00 . 
Furthermore, it is convenient to group the basis states into manifolds of similar energy, 
that we call n*-manifolds. Each consists of the three quantum defect states | n* + 3, 0, 0), 

| n* + 2,1, 0), | n* + 1, 2, 0) and the n* — 3 energetically degenerate states | n*, l > 2,0). 
For practical purposes and without loss of generality, we can neglect states with m^O. 
The projection of the electronic angular momentum onto the internuclear R axis parallel 
to the 2 -axis is a good quantum number. States with m^O have vanishing density on 
the ^-axis and hence show no interaction via the potential V. 

The resulting BO potential energy curves are shown in Fig. [D Using only the 
n* = 32 manifold gives the BO energy E 0 (cyan dots), including the n*-manifolds with 
\n* — 321 < 1 leads to the BO energy E\ (light blue dots) and taking \n* — 321 <5 results 
in the BO energy E 5 (dark blue dashed-dotted line). Here every increase in the size of 
the basis set leads to considerably lower potential curves. Surprisingly, even for large 
basis sets consisting of up to 10 n*-manifolds, the potential curves do not converge. 

The following calculation shows that this non-convergence is not a numerical 
artifact but inherent to the bare delta function potential used for the electron-perturber 
interaction and the diagonalization procedure. We will illuminate this in the following. 
Let us consider the finite subspace W spanned by N eigenstates | (pf) and investigate the 
diagonalization of H inside W in more detail. 

The projection of the Schrodinger equation 

HV(r) = EV( r) (10) 

onto W reads 

N 

E [*?*« + 2xa[%*(R)^(R)] a.j = Eoti 


( 11 ) 
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where we made the ansatz 

N 

^( r ) = • 

3 = 1 

Consequently, (1T21) will lead to the exact solution of (TTOl) in the limit N —y oo. 
From (HIT) we deduce 


( 12 ) 


(Xn 


|i^r =W V ‘ (R 


3 = 1 


E? - E 


(13) 


where J\f is a normalization constant. Reinserting this into (fTT!) yields the characteristic 
equation 


N 


1 + 2na[k\ 


i=l 


\^m 2 

E? -E 


= 0 . 


(14) 


The roots E of (fl4l) are now exactly the eigenvalues that one would obtain by a numerical 
diagonalization of (1111) . However, in contrast to numerical diagonalization methods there 
is no need to perform the approximation (|6]) for the involved /c-vector in numerical root 
finding algorithms. Hence, (THl) is even more general. 

To see the effects of diagonalization inside large basis sets, we take the limit N —> oo. 
The sum in (1T41) can then be replaced by the Green’s function of H° 

OO 

>*( r 'M( r ) 


G°(r,r',E) = 


E? - E 


(15) 


i= 1 


evaluated in the limit r = r' = R, where the sum above includes discrete and continuum 
states. However, for r —> r' the Green’s function diverges as 

G°(r,r',E) oc ^ ■ (16) 

Therefore the characteristic equation dlTl) possesses for IV — > oo no well defined solution 
E. This explains the nonconvergent behavior that we encounter with increasing size of 
the basis set for the diagonalization, as illustrated in Fig. [D 


4. Regularized delta function potential 


In the literature it is well known that the usefulness of a delta function potential in three 
dimensions is restricted since its action on irregular wave functions is not well defined 
and its spectrum is unbounded from below [251126] , We employ therefore the regularized 
delta function potential 

d 

R(r,R) = 27ra[k(R, E)}5{r ~ R)— p, (17) 

whose solutions fulfill the condition (|2]) exactly [27] , Using again the basis functions ipi 
we obtain 


0 = (£» 


E)oti + 2 r Ka[k(R 1 E)\tp* (R) 



(18) 
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Similar to (fTTl) and as shown in |2H] we derive that a nontrivial solution of (TTH1) exists 
only for energies satisfying 


1 + 2Tra[k(R, E )] 



7h*(R)A 

E? - 



= 0 


(19) 


r >R 


In contrast to (0711) this equation is not restricted to finite regions of the Hilbert space 
and yet well defined. It is crucial to carry out the summation before taking the derivative 
in (fTHD and ([TUI) as the opposite, i.e. shifting the derivative into the sum, would again 
lead to irregularities, due to (fT6l) . To this end, we insert explicitly the Green’s function 
(fT5l) into (TT9l) which yields 


1 + 2Tra[k(R, if)]G) ) eg (R, R, E) — 0 , 


( 20 ) 


where 


G?« g (R. R, E) 


' d_ 
dp 


(pG°( r,R,B)) 


r->n 


( 21 ) 


This is exactly the result obtained alternatively via the quantum defect theory in | !20lf29] 
and agrees nicely with the expectation that the pseudopotential approach and the 
quantum defect approach should result in equivalent potential curves. Eq. (1201) can be 
evaluated numerically by expressing the Green’s function in (1T51) in terms of Whittaker 
functions, which is possible for arbitrary quantum defects [30]. The BO potential energy 
curve E reg (R) obtained in this manner for the Rb(35s)-\-Rb(5s) state with the parameters 
(J9]) is shown (orange dashed curve) in FigUl Compared to the first order perturbation 
theory results, the outer minimum is around 30% deeper and shifted towards smaller R. 
This is an indication of the contribution of higher /-states and explains why the observed 
binding energies in pi] were larger than expected from the first order perturbation 
approximation (JH]). In [2], this was compensated partially by the introduction of an 
effective scattering length a e g[0]. 

More analytical insight is gained by performing first order perturbation theory with 
urn Additional care is necessary when using the regularized delta function potential 
in standard perturbation theory, due to the non-commutativity of summation and 
differential operations m Nevertheless, the first order corrections to the energy 
E® can be obtained in the normal manner as 


e\ 1] ’ = 2na[k{R, E°] j d 3 r<p*(r)5(r - R)^p^(r) 
= 2na[k(R, E°] |^(R)| 2 


( 22 ) 


where energetically degenerate (pi are chosen such that they diagonalize V in the 
degenerate manifold. Hence the first order energies of the regularized delta function 
potential (1UU1) coincide with the first order energies of the bare delta function potential 
(JH]). Therefore, even though we showed that the bare delta function potential is in 
general not adequate, it serves as a valid first order approximation for the regularized 
delta function potential. 
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More generally, Fig. |T| suggests that the bare delta function potential may even 
approximate the regularized delta function potential beyond the first order perturbation 
theory, by including only a few more, but not too many, adjacent states. For example 
the BO potential energy curve E 0 obtained by including only states adjacent to the 
n = 32 manifold, lies close to E reg (R). Here a more rigorous comparison between the 
gain in accuracy and the error accumulated when increasing the size of the basis set, 
would be desirable and should be carried out in the future. This would be especially 
useful as the regularized delta function potential can not be implemented in a numerical 
diagonalization procedure that relies on the basis functions <p; which are regular at R. 
We see this from (11811 . which, in contrast to dill) , can not be converted into a system 
of linear equations for the coefficients cq because the sum has to be carried out before 
taking the differential. 

5. Summary & Conclusion 

Aspects of the numerical implementation of the Fermi pseudopotential approach, 
relevant to the calculation of the BO potential energy curves of the ultralong-range 
Rydberg molecules bound by pure s-wave scattering are delineated. It is shown that 
the bare delta pseudopotential leads to unphysical results that diverge with increasing 
basis set size in any diagonalization scheme. This behavior is described analytically. 
The convergent BO potential curves, that agree exactly with the results obtained within 
the quantum defect theory can be produced by employing a regularized delta function 
potential. Although the bare delta function potential still yields the correct binding 
energies in the first order perturbation theory, the example of the Rb(35s) + Rb(5s) 
molecule indicates that beyond the first order effects are not negligible. By including 
some of the neighboring Rydberg basis states, some of these corrections can be 
reproduced. However it remains to be seen how many additional states can rigorously be 
included in exact diagonalization schemes and whether the zero-range pseudopotential 
for the p-wave-interaction p9], will require similar regularization as well. 
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